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Abstract 

One  challenge  associated  with  the  simulation  of  buried  detonations  involves  the  treatment 
of  the  multiphase  flow  phenomena  related  to  the  soil.  At  the  moment  when  the  soil  “shatters” 
into  a  dense  particulate  cloud  with  detonation  products  escaping  through  the  soil  particles,  the 
continuum  model  that  assumes  a  single  velocity  shared  by  the  blast  gas  and  the  soil  at  any  given 
point  no  longer  holds.  Instead  momentum  coupling  between  the  two  phases  has  to  be  modeled.  One 
characteristic  at  the  stage  of  soil  breaking  is  that  the  soil  fragments  packed  in  a  tight  configuration 
under  large  pressure  provide  significant  blockage  effect  characterized  by  large  particle  volume 
fractions.  Unfortunately,  traditional  drag  laws  do  not  address  the  momentum  coupling  between 
gas  and  solid  phase  under  the  condition  of  particle  high  volume  fraction  in  high  speed  blast  flows.  In 
order  to  develop  a  phenomenological  drag  model  to  characterize  the  momentum  coupling  between 
the  detonation  gas  and  soil  fragments  when  the  soil  initially  breaks  into  a  dense  particulate  cloud, 
we  conducted  a  series  of  numerical  simulations  on  the  scale  of  soil  fragments  by  only  considering 
a  small  region  occupied  by  a  mixture  of  blast  gas  and  soil  fragments  (so-called  particle-scale 
simulations).  A  drag  database  was  constructed  based  on  the  drag  force  collected  from  the  particle- 
scale  simulations  under  the  conditions  of  various  soil  volume  fractions  and  particle  sizes.  A  new 
drag  law  was  developed  using  data  regression  technique  to  characterize  the  dependency  of  the  drag 
force  exerted  on  particles  as  a  function  of  particle  volume  fraction  and  Reynolds  number  based  on 
particle  size.  The  proposed  drag  law  provides  satisfactory  representations  of  the  simulation  data, 
and  converge  to  traditional  drag  model  for  isolated  particles  when  the  particle  volume  fraction 
approaches  to  zero.  Finally,  we  applied  the  new  phenomenological  drag  law  in  the  simulations 
of  buried  detonation  at  the  stage  of  soil  breaking,  and  compared  the  results  with  the  ones  from 
continuum  model  in  which  the  gas  and  the  soil  are  assumed  moving  at  the  same  speed. 

key  word:  Drag  law,  Computational  fluid  dynamics,  Dense  particulate  flow,  Blast  wave,  Data 


REGRESSION 


UNCLASSIFIED:  Distribution  Statement  A.  Approved  for  public  release 


3 


1  Introduction 

In  CFD  simulations  of  explosive  blast  wave  traveling  through  broken  soil  fragments  in  the  case  of 
buried  detonations,  while  the  behavior  of  the  blast  wave  is  simulated  in  the  regular  flow  solver,  the 
small  scale  and  large  number  of  soil  fragments  make  it  impractical  to  simulate  the  behavior  of  soil 
particles  directly  in  the  flow  solver.  Therefore,  multiphase  coupling  effects  between  solid  phase  (soil 
fragments)  and  gas  phase  (detonation  gas)  need  to  be  modeled.  The  momentum  transfer  between  two 
phases  is  one  of  essential  coupling  terms  to  be  considered  since  it  is  critical  for  the  accurate  modeling  of 
blast  loading  on  the  target.  One  challenge  associated  with  momentum  transfer  between  two  phases  for 
buried  detonation  is  how  to  correctly  model  the  behavior  of  the  soil  that  surrounds  the  expanding  cloud 
of  detonation  products.  At  the  moment  when  soil  “shatters”,  or  breaks  into  a  cloud  of  particulates, 
the  continuum  model  that  assumes  a  single  velocity  shared  by  the  gas  and  soil  is  no  longer  appropriate. 
Instead,  a  drag  model  that  describes  the  functional  dependency  of  the  drag  force  exerted  on  soil  is 
needed.  One  characteristic  at  the  initial  stage  of  soil  breaking  is  that  the  soil  fragments  are  packed 
in  a  tight  configuration  under  large  pressure,  providing  a  significant  blockage  effect  characterized  by 
large  volume  fractions. 

There  have  been  significant  amount  of  effort  on  investigating  particle  drag  models  in  past  decades, 
ranging  from  theoretical  work  in  Stokes  flow  to  empirical  formulations  [Crowe,  1998].  Traditional 
multiphase  models  normally  consider  a  dispersed  phase  where  drag  laws  are  developed  assuming  that 
each  particle  is  isolated.  As  the  number  of  particles  increases  in  a  fixed  space,  two  effects  become 
important:  1)  particle  drag  laws  can  no  longer  assume  isolated  particles  as  one  particle  will  be  in  the 
viscous  wake  of  other  particles,  and  2)  a  non-negligible  amount  of  volume  occupying  by  the  particles 
will  cause  non-trivial  blockage  effects  in  the  fluid  phase.  Collectively,  these  two  effects  are  considered 
as  large  volume  fraction  effects,  and  they  only  occur  when  particle  volume  fraction  becomes  significant 
(greater  than  1%).  However,  there  is  little  information  available  on  the  drag  of  particles  in  dense  par¬ 
ticle  clouds,  with  mostly  commonly  mentioned  work  coming  from  the  sedimentation  and  fluidization 
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studies  [Gidaspow,  1994],  which  is  not  applicable  to  the  high  speed  blast  flows.  Properly  capturing 
large  volume  fraction  effects  is  required  to  reliably  describe  the  dense  particulate  flows  in  blast  wave 
that  occur  in  buried  detonations.  However,  current  available  models  make  no  attempt  to  capture  these 
effects.  The  goal  of  this  study  is  to  develop  appropriate  phenomenological  drag  laws  with  inclusion  of 
the  effects  of  large  volume  fractions  to  correctly  model  multiphase  effect  in  buried  detonation  simula¬ 
tions  at  the  stage  of  soil  failing.  In  this  study,  we  carried  out  the  numerical  simulations  in  the  scale  of 
soil  fragment  by  zooming  into  the  a  small  region  occupied  by  the  mixture  of  soil  fragments  and  the  gas, 
and  the  drag  force  acting  on  each  soil  fragments  was  computed  in  this  particle-scale  simulation.  The 
database  that  is  a  collection  of  the  drag  force  from  a  series  of  particle-scale  simulations  under  various 
particle  volume  fractions  and  sizes  was  established.  Then  the  drag  laws  were  obtained  by  analyzing 
the  data  from  the  database.  In  the  following,  the  first  step  of  particle-scale  numerical  simulations 
including  geometry  modeling  and  mesh  generation  will  be  described,  followed  by  the  presentation  of 
the  simulation  set-up  and  results.  The  database  of  collected  drag  force  and  the  drag  laws  based  on  the 
regression  of  the  database  will  be  discussed  in  the  end. 

2  Geometry  model  and  mesh  generation 

We  wish  to  simulate  the  coupling  effects  between  gas  and  particles  under  high  volume  fractions  in 
the  moments  associated  with  soil  failure.  The  geometry  of  the  particles  in  this  circumstance  cannot 
be  approximated  as  spherical  objects  due  to  the  high  packing.  In  order  to  obtain  very  large  volume 
fractions,  we  need  to  construct  particle  geometries  such  that  the  particles  can  be  expanded  and  returned 
to  their  original  continuous  state.  To  accomplish  this  task,  we  modeled  the  particles  by  constructing  a 
Voronoi  graph  about  a  random  set  of  seed  points.  The  Voronoi  graph  mimics  cleaves  that  form  between 
particles, and  expansion  of  the  the  Voronoi  graph  creates  gaps  where  the  passage  of  gas  was  modeled. 
The  Voronoi  diagram  subdivides  a  region  into  a  set  of  Voronoi  cells,  each  of  which  is  associated  with 
a  specific  “seed”  point.  A  Voronoi  cell  has  the  property  that  all  points  contained  within  the  cell  are 
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Figure  1:  Two-dimensional  represnetation  of  Voronoi  diagram  (image  from 

http:  / /www. donhavey.com/blog/ tutorials /tutorial-7-voronoi-diagrams/) . 


closer  to  its  seed  point  that  to  all  other  seed  points  in  the  domain.  Figure  1  shows  a  Voronoi  diagram 
for  a  planar  set  of  randomly  generated  points.  The  “borders”  of  each  cell  have  equal  distance  between 
the  seed  points  of  adjacent  cells.  This  creates  convex  regions  that  can  be  organized  or  grouped  as  soil 
fragments.  In  order  to  reduce  the  size  of  the  computational  domain  needed  to  simulate,  the  Voronoi 
graph  was  constructed  to  be  periodic.  This  approach  accommodates  creating  geometry  of  arbitrarily 
large  particle  volume  fractions. 

Four  pieces  of  information  were  used  to  create  the  simulated  particle  formation,  namely  particle 
centers  (random  points),  Voronoi  points,  Voronoi  neighbors  around  the  Voronoi  points,  and  the  Voronoi 
points  for  each  Voronoi  region.  The  particles  were  assembled  by  isolating  the  neighbor  map  for  each 
Voronoi  region  and  only  the  regions  with  positive  neighbors  and  overlapped  with  the  unit  cubes  were 
considered.  The  cube  and  the  particles  were  then  displaced/expanded  from  the  center  of  the  unit 


UNCLASSIFIED:  Distribution  Statement  A.  Approved  for  public  release 


6 


cube  like  expanding  the  universe  to  create  various  volume  fraction.  The  particle  faces  were  then 
formed  by  traversing  the  particle  neighbor  map  and  only  considering  points  that  lie  within  a  plane  to 
form  a  closed  loop  that  represents  the  surface  of  the  particle.  The  faces  were  then  intersected  with 
Y  and  Z  planes  of  the  cubes.  Any  parts  that  were  not  fully  enclosed  in  the  cube  were  trimmed. 
The  functionality  of  the  SolidMesh  library  [Gaither  et  al.,  2000]  was  used  to  assemble  the  geometric 
description  and  perform  any  necessary  geometric  operations.  The  resultant  geometry  was  assembled 
into  a  solid  modeling  representation  and  output  as  a  SolidMesh  save  state  which  can  be  read  into 
SolidMesh  for  mesh  generation. 

There  are  competing  concerns  that  must  be  balanced  in  building  geometry  model.  We  would  like 
to  have  as  many  soil  fragments  as  possible  to  better  approximate  an  actual  particle  cloud;  however,  we 
cannot  afford  to  have  meshes  that  contain  a  large  number  of  cells  due  to  computational  cost.  Based  on 
some  preliminary  studies,  we  determined  that  having  10  soil  fragments  per  unit  volume  (with  overall 
number  of  fragments  at  about  40)  and  utilizing  periodic  boundaries  provide  a  good  compromise  between 
computational  cost  and  the  ability  to  capture  true  multiphase  effects.  Figure  2  shows  the  fluid  portion 
of  the  periodic  cube  of  particulate  grain.  In  this  figure,  the  outer  faces  on  X  direction  are  placed  a  few 
unit  away  from  the  expanding  unit  cube  to  decrease  the  boundary  effect  for  the  flow  simulations.  In 
this  model,  10  random  points  in  one  unit  volume  are  generated  (with  total  of  34  random  points)  and 
expansion  rate  is  1.25  which  corresponds  to  particle  volume  fraction  of  51.2%.  Only  those  particles 
that  are  fully  enclosed  in  the  cube  and  that  intersected  the  Y  and  Z  planes  were  considered  in  the 
simulations.  The  particles  with  surrounding  faces  turned  off  are  shown  in  figure  3. 

3  Numerical  Simulations 

When  the  flow  is  dilute  where  the  particle  volume  fraction  is  small,  the  drag  coefficient  exerted  on 
the  isolated  particle  is  a  function  of  Reynolds  number  based  on  particle  size  and  the  relative  velocity 
between  two  phases,  according  to  traditional  drag  laws  [Clift  and  Gauvin,  1970].  Thus  the  particle 
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Figure  2:  Outer  face  of  periodic  cube  of  particle  cloud 


Figure  3:  The  underlying  particle  cloud 
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size  plays  a  dominant  role  in  the  cases  of  low  particle  volume  fraction  where  the  viscous  effects  are 
important.  When  the  particles  are  densely  packed  in  the  flow  field,  the  drag  between  particles  and  the 
gas  mainly  comes  from  inviscid  nature  of  the  flow,  and  the  particle  volume  fraction  plays  a  significant 
role.  In  order  to  evaluate  the  importance  of  soil  size  and  soil  volume  fraction  on  momentum  transfer 
between  soil  fragments  and  high  pressure  detonation  gas,  which  leads  to  a  phenomenological  drag  model 
of  blast  gases  flowing  through  particles,  a  series  of  cases  with  high  pressure  flows  passing  through  a 
number  of  soil  fragment  pieces  with  different  particle  sizes  and  volume  fractions  was  simulated  in  the 
small  zoomed-in  region  of  the  domain  (so-called  particle-scale  simulations).  About  40  pieces  of  soil 
fragments  existed  in  the  simulation  domain.  There  were  six  different  expansion  ratios  of  soil  fragment 
centroids,  namely  1.2,  1.25,  1.5,  2,  5  and  10,  which  corresponded  to  soil  volume  fractions  of  57.8%, 
51.2%,  29.6%,  12.5%,  8%  and  0.1%.  This  range  of  volume  fractions  covers  the  cases  with  densely  packed 
soil  blocks  to  the  situations  with  dilutely  distributed  soils.  Figures  4  and  5  show  the  configurations 
of  soil  fragments  with  highest  and  lowest  particle  volume  fractions  in  our  particle-scale  simulations. 
Under  each  volume  fraction,  five  different  soil  sizes  are  evaluated,  which  is  on  the  order  of  0.1mm, 
0.3 mm,  1  mm,  3 mm  and  1cm.  The  different  soil  size  is  achieved  by  setting  corresponding  reference 
length  in  the  simulations.  The  computational  domain  is  a  rectangular  box  with  periodic  boundary 
conditions  imposed  at  four  boundaries  surrounding  the  soil  fragments  in  both  Y  and  Z  directions 
simulating  an  infinite  field  of  debris.  The  flow  direction  is  along  the  X  axis  and  we  provided  a  buffer 
region  ahead  of  and  behind  the  flow  so  that  the  flow  could  equilibrate  before  reaching  the  boundaries. 

In  the  case  of  blast  simulations,  the  gas  flow  and  particles  are  in  complete  momentum  equilibrium 
and  hence  have  the  same  velocity  at  the  moment  when  particle  cloud  is  formed.  The  gas  is  able  to 
pass  through  gaps  between  the  soil  fragments.  Note  that  while  the  gas  and  particles  may  be  traveling 
at  supersonic  speed,  the  relative  velocities  between  the  two  materials  is  subsonic  at  the  moment  of 
soil  failure.  Thus  the  numerical  simulations  focus  on  flow  through  the  particulate  cloud  using  subsonic 
boundary  conditions.  Also,  since  at  the  time  of  particle  formation  the  gas  density  is  typically  much 
smaller  than  the  solid  material,  the  velocity  gained  by  the  particles  is  very  small  relative  to  the  gas 
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Figure  4:  Configuration  of  soil  fragments  with  highest  particle  volume  fraction  (57.8%) 


Figure  5:  Configuration  of  soil  fragments  with  lowest  particle  volume  fraction  (0.1%) 
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velocity.  Therefore,  it  is  reasonable  to  assume,  for  the  purpose  of  developing  drag  laws,  that  the 
particles  are  stationary.  While  the  particles  will  move,  their  movements  will  be  over  a  much  larger 
time  scale  than  the  time  scale  of  the  gas  flow.  A  high  temperature  gas  flows  through  soil  fragments. 
The  current  model  employed  a  gas  modeled  as  an  ideal  gas  with  temperature  of  2000 A'.  The  inflow 
boundary  condition  was  isentropic  with  total  pressure  of  10 atm.  There  were  several  number  of  pressure 
set  up  at  outflow  boundary,  namely  2 atm,  3atm,  latm,  5atm  and  7 atm.  Simulations  from  these 
different  outflow  pressure  suggested  that  flow  undergoes  choking  due  to  the  blockage  effect  of  soils,  as 
shown  in  figure  6.  Therefore,  drags  exerted  on  particles  were  not  sensitive  to  outflow  pressure  except 
for  those  closest  to  the  outflow  boundary.  This  conclusion  was  also  justified  in  the  simulation  results 
under  different  outflow  pressures.  Thus,  we  only  chose  the  simulations  with  outflow  boundary  pressure 
of  5 atm  for  construction  of  our  drag  database.  The  pressure  distributions  on  the  soil  fragments  under 
different  particle  volume  fraction  are  demonstrated  in  figs  7  to  12.  It  is  clear  that  the  pressure  change 
across  the  domain  occupied  by  soil  blocks  decreases  as  the  the  flow  becomes  more  dilute. 

Two  non-dimensional  variables  needed  in  the  drag  model  are  drag  coefficient  {Cd),  Reynolds  number 
based  on  particle  size  ( Red ).  They  are  defined  as  : 


Cd 


Red  = 


_  2  *  Fd 
Au 2 

p(u  —  up)d 


(1) 

(2) 


where  Fd  is  the  drag  force  exerted  on  soil  fragments,  while  u  and  up  denote  the  speed  of  gas  and  soil. 
As  mentioned  before,  up  is  considered  zero  in  our  simulations.  The  density  and  dynamic  viscosity  of 
the  gas  are  represented  by  p  and  p,  respectively.  The  variable  A  represents  the  projected  area  of  particle 
in  the  flow  direction,  while  d  represents  the  diameter  of  the  particle.  In  our  simulations,  the  shape  of 
soil  fragments  was  non-spherical,  thus  particle  diameter,  d,  is  defined  by  the  volume-equivalent-sphere 
diameter  or  nominal  diameter  represented  as  d  =  sJUVJtx  where  V  is  the  particle  volume. 


UNCLASSIFIED:  Distribution  Statement  A.  Approved  for  public  release 


11 


Figure  6:  Mach  number  distribution  at  z  =  0  cutting  plane 


1.000e+06  r 

8.250e+05  I 
6.500e+05  I 
4.750c+05  I 
3.000c+05  8 


Figure  7:  Pressure  distribution  on  soil  fragments  for  57.8%  particle  volume  fraction 
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Figure  8:  Pressure  distribution  on  soil  fragments  for  51.2%  particle  volume  fraction 
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Figure  9:  Pressure  distribution  on  soil  fragments  for  29.6%  particle  volume  fraction 


Figure  10:  Pressure  distribution  on  soil  fragments  for  12.5%  particle  volume  fraction 
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Figure  11:  Pressure  distribution  on  soil  fragments  for  8%  particle  volume  fraction 


Figure  12:  Pressure  distribution  on  soil  fragments  for  0.1%  particle  volume  fraction 
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Since  there  is  a  wide  distribution  of  values  of  the  variables  in  the  simulation  domain,  it  is  important 
to  define  a  sensible  set  of  variables  for  nondimensionalization.  Since  in  practice  we  will  be  utilizing  this 
database  to  represent  the  integrated  effects  throughout  a  cell  volume  that  is  generally  much  larger  than 
soil  particulate  sizes,  we  adopted  values  that  result  from  integrations  around  the  region  of  each  particle 
in  order  to  obtain  values,  such  as  fluid  velocity  and  density,  that  are  used  for  nondimensionalization. 
To  perform  this  integration  we  implemented  a  tool  that  computes  the  volume-averaged  local  variables 
(e.g.  u,p,p)  in  the  median  volume  around  each  soil  fragment,  defined  as: 


J2  SjVoli 

J2  v°h 


(3) 


where  s  stands  for  any  median  volume  averaged  variable,  and  vol  is  the  volume  of  the  cell.  Index  i 
represents  the  cell  that  within  the  median  volume  of  the  soil  fragment.  The  median  volume  associated 
with  soil  fragment  is  all  the  cells  that  are  closer  to  the  given  soil  fragment  than  any  others.  It  was 
observed  that  there  were  a  significant  amount  of  noises  in  the  data  around  each  soil  fragment  due  to 
localized  variability  in  the  particle  arrangement  and  shape.  Therefor  we  determined  to  use  averaged 
values  over  the  computational  domain  to  compute  drag  coefficients  for  the  database.  The  averaged 
values  of  variables,  such  as  drag  force,  velocity,  diameter  and  other  properties  of  soil,  were  obtained  by 
averaging  the  corresponding  information  associated  with  each  soil  fragment  excluding  the  parts  whose 
median  volumes  extended  all  the  way  to  the  inflow  or  outflow  boundaries,  i.e.  particles  that  were  at 
the  leading  or  trailing  edges  of  the  simulated  particles. 

Note  that  the  drag  force  on  soil  fragment  in  the  case  of  high  particle  volume  fraction  comes  from  the 
combination  of  the  local  effect  of  flow  around  the  soil  fragments  and  the  global  effect  associated  with 
the  bulk  pressure  change  due  to  the  blockage  presented  by  the  large  volume  fraction.  When  particle 
volume  fraction  decreases,  the  pressure  change  across  the  computational  domain  decreases,  thus  the 
drag  force  becomes  dominated  by  local  inviscid  and  viscous  effects.  This  distinction  is  important  since 
in  the  context  of  a  drag  law,  the  mesh  in  the  nominal  CFD  solution  is  much  larger  than  the  particle 
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size  so  the  bulk  effects  of  pressure  gradient  acting  over  the  surface  of  the  particles  will  be  directly 
captured  as  an  average  behavior  by  the  CFD  solver,  whereas  local  effects  need  to  be  provided  by  a 
phenomenological  model.  To  obtain  this  distinction  between  average  bulk  effects  resolved  by  the  large 
scales  and  the  small  scale  effects,  we  consider  the  effect  of  the  bulk  pressure  gradient  in  these  detailed 
particle-scale  simulations.  In  order  to  characterize  this  local  effect,  we  defined  another  drag  coefficient, 
Ccii ,  that  excludes  the  global  pressure  gradient,  as  follows 


Cdi 


2  *{Fd~  Fpg) 
Au 2 


(4) 


where  Fpg  is  the  force  due  to  the  bulk  pressure  gradient.  This  bulk  pressure  gradient  is  computed  by 
integrating  pressure  over  the  surface  of  median  volume  which  excludes  the  particle  local  effect  such  as 
stagnation,  and  then  by  scaling  the  value  to  the  area  of  corresponding  soil  fragment  based  on  the  ratio 
of  the  surface  area  of  the  soil  fragment  to  the  surface  area  of  the  median  volume. 


4  Database  and  Drag  Model 

The  database  of  drag  is  presented  in  table  1.  As  mentioned  previously,  the  values  in  the  database 
are  averaged  over  computational  domain  due  to  the  large  fluctuations  in  the  data  around  each  soil 
fragment.  The  magnitude  of  Reynolds  number  covers  from  102  to  105,  and  particle  volume  fraction 
has  a  2  order  range  of  values.  In  the  table,  CdP  is  the  parameter  based  on  the  bulk  pressure  effect  on 
particles,  and  defined  as  Cdp  =  Cd  —  Cdi •  It  shows  that  the  contribution  from  bulk  pressure  gradient 
decreases  as  flow  becomes  more  dilute  due  to  less  blockage  effect  from  the  expanding  soil  fragments. 
When  the  particle  volume  fraction  approaches  zero,  the  effect  from  bulk  pressure  gradient  imposed  on 
soils  is  negligible,  and  the  drag  force  results  from  the  local  effect  around  soils. 

We  applied  data  regression  to  obtain  local  drag  coefficient  dependency  on  particle  volume  fraction 
and  Reynolds  number  based  on  particle  size,  i.e.  Cdi  =  f(Red,oi).  It  is  observed  in  table  1  that  drag 
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coefficient  based  on  the  total  drag  force  ( Cd )  increases  monotonically  with  particle  volume  fraction, 
while  local  drag  coefficient  that  excludes  the  bulk  pressure  term,  Cdi ,  does  not  have  monotonic  behavior 
with  particle  volume  fraction.  However,  the  non-dimensional  parameter  that  characterizes  the  bulk 
pressure  gradient  ( CdP  =  Cd  —  Cdi)  demonstrates  monotonic  increase  with  particle  volume  fraction. 
Due  to  the  lack  of  the  monotone  behavior  of  Cdi ,  we  determined  that  a  better  strategy  was  to  perform 
data  regression  on  Cd  and  CdP  such  that  a  more  accurate  fitting  would  be  achieved.  The  model  for  the 
local  drag  coefficient  (Cdi)  was  then  computed  by  subtracting  contribution  due  to  the  bulk  pressure 
gradient  (CdP)  from  the  total  drag  coefficient  (Cd). 

A  data  regression  open-source  software  package,  R-language  (http://www.r-project.org/),  which 
has  been  widely  used  in  statistical  community  was  adopted  as  our  data  analysis  toolset.  The  new 
drag  model  should  converge  to  a  traditional  drag  law  for  isolated  particles  as  particle  volume  fraction 
approaches  to  zero.  One  of  the  most  widely  used  drag  law  for  isolated  spheres  is  the  empirical  equation 
of  Clift  and  Gauvin.  [Clift  and  Gauvin,  1970],  which  is  valid  for  particle  Reynolds  number  up  to  3  x  105 
and  reads  as: 


Cd0  =  —  (1  +  0.15i?e°'687)  + - — - ^ 

Red V  d  1  1  +  4.25  x  104i?e7116 


(5) 


where  Cdo  stands  for  the  drag  coefficient  when  the  particle  volume  fraction  is  zero.  In  this  study, 
two  additional  effects  needs  to  be  taken  into  account  in  the  drag  law  of  isolated  particles,  namely 
compressibility  effect  and  non-spherical  particle  effect  .  When  the  particle  Reynolds  number  is  not 
too  low  (Red  >  45),  the  formation  of  shock  waves  on  the  particles  and  attendant  wave  drag  increase 
the  drag  force  acting  on  the  particles.  This  compressibility  effect  becomes  significant  when  Mach 
number  reaches  0.6,  which  is  the  critical  Mach  number.  We  adopted  the  drag  model  developed  by 
Loth  [Loth,  2008]  for  the  correction  of  the  Mach  number  effect: 


Cd0 


—  (1  +  0.15i?e°/87)  *  (1-0 

tied 


0.258  *  CM 
1.0  +  514*  GAT 


0.42  *  CM 

1  +  4.25  x  104  *  CM  *  Re~1Ae 


(6) 
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Table  1:  Drag  coefficients  vs  Reynolds  number  (Red)  and  volume  fraction  (a)  under  different  particle 
sizes  _ 


cd 

Cdi 

cd  -Cdl 

Red 

a 

7.798 

3.239 

4.559 

148.9 

0.579 

4.600 

2.190 

2.410 

188.9 

0.512 

0.1  mm 

1.950 

1.298 

0.652 

333.0 

0.296 

1.454 

1.195 

0.259 

449.9 

0.125 

0.973 

0.950 

0.023 

719.7 

0.008 

0.909 

0.903 

0.006 

834.0 

0.001 

2.930 

1.167 

1.763 

710.3 

0.579 

2.284 

1.041 

1.243 

787.6 

0.512 

0.3  mm 

1.532 

0.992 

0.540 

1125.4 

0.296 

1.317 

1.070 

0.247 

1373.0 

0.125 

0.942 

0.919 

0.023 

2161.0 

0.008 

0.853 

0.847 

0.006 

2502.1 

0.001 

1.897 

0.712 

1.185 

2883.6 

0.579 

1.685 

0.745 

0.940 

3028.3 

0.512 

dr 

-r  1mm 

1.440 

0.920 

0.520 

3830.2 

0.296 

1.293 

1.045 

0.248 

4578.5 

0.125 

0.953 

0.930 

0.023 

7154.3 

0.008 

0.860 

0.854 

0.006 

8346.7 

0.001 

1.685 

0.613 

1.072 

9104.4 

0.579 

1.515 

0.660 

0.855 

9576.1 

0.512 

dr 

-r  3 mm 

1.391 

0.869 

0.522 

11602.8 

0.296 

1.300 

1.051 

0.249 

13547.7 

0.125 

0.938 

0.916 

0.022 

21514.5 

0.008 

0.822 

0.816 

0.006 

24971.0 

0.001 

1.567 

0.560 

1.007 

31335.5 

0.579 

1.511 

0.691 

0.820 

33146.0 

0.512 

d ■ 

~  1  cm 

1.369 

0.853 

0.516 

38640.9 

0.296 

1.320 

1.072 

0.248 

44995.2 

0.125 

0.974 

0.951 

0.023 

71603.3 

0.008 

0.859 

0.854 

0.005 

83425.2 

0.001 
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where  CM  and  CM  are  function  of  Mach  number,  M,  and  are  defined  as: 


CM  =  7T  *  tanh(3.0  *  log(M  +  0.1))  ifM  <  1.45  (7) 

O 

CM  =  2.044  +  0.2  *  exp(— 1.8  *  (log{^-))2)  ifM  >  1.45  (8) 

GM  =  1.0  —  1.525  *  M4  ifM  <0.89  (9) 

GM  =  0.0002  +  0.0008  *  tan/i(12.77  *  (M  -  2.02))  ifM  >0.89  (10) 


For  non-spherical  particles  in  an  infinite  medium,  it  is  generally  recognized  that  a  shape  factor  must  be 
included  in  the  expression  of  drag  coefficient.  Shape  factor  of  Tran-Cong  et.  al  [Tran-Cong  et  ah,  2004] 
was  utilized  to  modify  the  drag  model  as  follows: 


Cda  =  Cdo  *  (— r)2  (11) 

a 

where  Cda  represents  the  drag  coefficient  with  consideration  of  particle  irregularity,  and  d  is  the  volume- 
equivalent-sphere  diameter  defined  earlier.  The  surface-equivalent-sphere  diameter,  dj\ ,  is  written 
as  6?4  =  ^jAAp/ir  where  Ap  is  the  projected  area  of  the  particle.  The  ratio  is  used  as  the 

dimensionless  shape  factor  to  account  for  the  non-spherical  effect  of  particles  on  drag  force. 

Using  muti-variable  nonlinear  data  regression  technique  in  R-language,  the  new  drag  models  with 
particle  volume  fraction  being  a  dependency  variable  were  obtained: 

Cd  =  Cda  +  Al  *  a  *  (1.0  +  j^^ci  )4  *  ot  *  (1.0  +  p^ci )  (12) 

Cdi  =  Cda+AUa*(lA+^)4+BUa*(l.0+^)-A2*a*(1.0+-^)6-B2*a*(l.0+-^)  (13) 

Red  Red  Red 
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where  fitting  coefficients  are  given  by  A1  =  0.5048,51  =  0.9858,  cl  =  0.5707,51  =  34.8  for  the  total 
drag  coefficient,  Cd,  and  A2  =  1.876,52  =  1.4285,  c2  =  0.4775,52  =  9.4926  for  the  pressure  drag 
coefficient,  CdP ■  It  is  obvious  that  the  drag  models  become  conventional  drag  law  of  isolated  particles 
when  particle  volume  fraction  is  zero.  The  drag  models  show  that  drag  force  increases  with  particle 
volume  fraction, and  the  steepest  increase  happens  at  small  a  and  Red,  which  was  evaluated  by  the 
higher  power  in  the  model  functions.  Figures  13  and  14  demonstrate  comparison  between  numerical 
data  and  the  results  from  our  new  drag  models  based  on  total  drag  and  local  drag  force,  respectively. 
In  the  figures,  the  surfaces  defined  by  the  green  grids  represent  the  data  from  particle-scale  simulations 
while  the  surfaces  from  the  drag  models  are  shown  in  blue.  Tables  2  and  3  show  comparison  of  drag 
coefficients  from  database  and  proposed  model.  The  relative  error  is  defined  as  the  absolute  value  of 
the  difference  between  numerical  data  and  result  from  the  model  divided  by  numerical  data.  It  is  seen 
that  the  fitting  errors  at  low  end  of  particle  volume  fraction  are  relatively  higher  with  our  models  due 
to  an  over-prediction  compared  with  the  database.  The  square  root  of  sum  of  square  of  each  relative 
error  is  9.4%  for  Cd,  and  10.2%  for  Cdi ■  These  results  indicate  that  the  proposed  drag  models  provide 
satisfactory  agreement  with  data  from  particle-scale  simulations. 

Before  we  close  this  section,  it  is  helpful  to  discuss  some  issues  not  concerned  in  the  development 
of  our  drag  model  so  that  we  understand  some  of  the  limitations.  In  our  particle-scale  simulations  in 
which  soil  clouds  uniformly  expanding  in  the  domain,  it  is  assumed  that  soil  fragments  remains  their 
orientations  during  the  simulations  while  in  reality  the  soil  pieces  might  rotate  due  to  the  moment 
imposed  on  them.  The  ignorance  of  the  particle  rotation  could  cause  the  inaccurate  prediction  of  the 
drag  force.  However,  the  complexities  involved  in  rotation  and  spinning  of  the  soil  fragments  would 
make  the  development  of  drag  model  intractable  in  the  current  stage,  and  our  present  drag  model 
mainly  reflect  the  blockage  effect  on  multiphase  flows.  Since  the  time  scale  of  blast  detonation  is  very 
small,  it  is  reasonable  to  treat  soil  fragments  irrotational.  Another  factor  not  considered  in  the  model 
is  compressibility  effect  under  high  particle  volume  fraction.  We  assume  that  flow  is  choked  and  Mach 
number  does  not  change  significantly  over  the  simulation  domain.  In  the  future  study,  database  will 


UNCLASSIFIED:  Distribution  Statement  A.  Approved  for  public  release 


20 


100 


100000 


Figure  13:  Comparison  of  variation  of  total  drag  coefficient  ( Cd )  with  particle  volume  fraction  and 
Reynolds  number  between  numerical  data  and  proposed  model 


Figure  14:  Comparison  of  variation  of  local  drag  coefficient  (C^i)  with  particle  volume  fraction  and 
Reynolds  number  between  numerical  data  and  proposed  model 
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Table  2:  Comparison  of  total  drag  coefficients  (Cd)  between  database  and  model  results 


Cd  from  database 

Cd  from  the  model 

relative  error 

volume  fraction 

7.798 

7.753 

0.6  % 

0.579 

4.600 

4.729 

2.8  % 

0.512 

1.950 

1.986 

1.9  % 

0.296 

1.454 

1.427 

1.9  % 

0.125 

0.973 

1.121 

15.2  % 

0.008 

0.909 

1.092 

20.1  % 

0.001 

2.930 

2.769 

5.5  % 

0.579 

2.284 

2.335 

2.2  % 

0.512 

1.532 

1.564 

2.1  % 

0.296 

1.317 

1.237 

6.1  % 

0.125 

0.942 

1.023 

8.6  % 

0.008 

0.853 

1.004 

17.8  % 

0.001 

1.897 

1.975 

4.1  % 

0.579 

1.685 

1.796 

6.6  % 

0.512 

1.440 

1.375 

4.5  % 

0.296 

1.293 

1.131 

12.5  % 

0.125 

0.953 

0.964 

1.2  % 

0.008 

0.860 

0.950 

10.5  % 

0.001 

1.685 

1.739 

3.2  % 

0.579 

1.515 

1.611 

6.4  % 

0.512 

1.391 

1.285 

7.6  % 

0.296 

1.300 

1.075 

17.3  % 

0.125 

0.938 

0.930 

0.9  % 

0.008 

0.822 

0.919 

11.7  % 

0.001 

1.567 

1.618 

3.3  % 

0.579 

1.511 

1.512 

0.0  % 

0.512 

1.369 

1.230 

10.2  % 

0.296 

1.320 

1.038 

21.4  % 

0.125 

0.974 

0.906 

7.1  % 

0.008 

0.859 

0.896 

4.3  % 

0.001 
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Table  3:  Comparison  of  local  drag  coefficients  (Cdi)  between  database  and  model  results 


Cdi  from  database 

Cdi  from  the  model 

relative  error 

volume  fraction 

3.239 

3.189 

1.5% 

0.579 

2.190 

2.364 

8.0% 

0.512 

1.298 

1.292 

0.5% 

0.296 

1.195 

1.157 

3.2% 

0.125 

0.950 

1.104 

16.2% 

0.008 

0.903 

1.090 

20.6% 

0.001 

1.167 

1.040 

10.9% 

0.579 

1.041 

1.069 

2.7% 

0.512 

0.992 

0.994 

0.2% 

0.296 

1.070 

1.005 

6.1% 

0.125 

0.919 

1.009 

9.7% 

0.008 

0.847 

1.003 

18.3% 

0.001 

0.712 

0.750 

5.4% 

0.579 

0.745 

0.810 

8.7% 

0.512 

0.920 

0.871 

5.4% 

0.296 

1.045 

0.922 

11.8% 

0.125 

0.930 

0.951 

2.3% 

0.008 

0.854 

0.948 

11.0% 

0.001 

0.613 

0.671 

9.5% 

0.579 

0.660 

0.726 

10.0% 

0.512 

0.869 

0.814 

6.3% 

0.296 

1.051 

0.879 

16.3% 

0.125 

0.916 

0.917 

0.2% 

0.008 

0.816 

0.917 

12.3% 

0.001 

0.560 

0.631 

12.6% 

0.579 

0.691 

0.682 

1.3% 

0.512 

0.853 

0.779 

8.7% 

0.296 

1.072 

0.849 

20.7% 

0.125 

0.951 

0.894 

6.0% 

0.008 

0.854 

0.894 

4.8% 

0.001 
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include  Mach  number,  and  drag  coefficient  on  soil  fragments  will  be  a  function  of  Mach  number  besides 
two  variables  in  the  present  work. 

5  Summary  and  concluding  remarks 

The  objective  of  this  study  was  to  develop  a  phenomenological  drag  model  to  characterize  the  momen¬ 
tum  coupling  between  the  detonation  gas  and  soil  fragments  when  the  soil  initially  breaks  into  a  dense 
particulate  cloud.  We  conducted  a  series  of  numerical  simulations  on  the  scale  of  soil  fragment  by 
only  considering  a  small  region  occupied  by  a  mixture  of  blast  gas  and  soil  fragments.  Voronoi  graph 
was  adopted  to  model  the  geometry  of  soil  breaking.  Different  expansion  ratios  between  the  Voronoi 
points  were  applied  to  form  six  soil  volume  fractions,  covering  the  cases  with  densely  packed  soil  cloud 
to  the  situations  with  dilutely  distributed  soil  fragments.  Five  different  particle  sizes  were  considered 
under  each  soil  volume  fraction.  In  the  simulations,  only  one  value  of  pressure  was  imposed  at  the 
outflow  boundary  since  the  flow  underwent  choking  due  to  the  blockage  of  the  soil.  A  drag  database 
was  constructed  based  on  the  drag  force  collected  from  the  particle-scale  simulations.  A  new  drag 
law  was  developed  using  data  regression  technique  to  characterize  the  dependency  of  the  drag  force 
exerted  on  particles  as  a  function  of  particle  volume  fraction  and  Reynolds  number  based  on  particle 
size.  The  irregularity  of  the  particles  was  also  considered  in  the  drag  model.  The  proposed  drag  law 
provides  satisfactory  representations  of  the  simulation  data,  and  converge  to  traditional  drag  model 
for  isolated  particles  when  the  particle  volume  fraction  approaches  to  zero. 
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